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ABSTRACT 


The  reliability  of  a  structural  system  at  a  particular  time 
depends  on  the  damage  level  in  the  system.  When  the  damage  level 
exceeds  a  critical  value,  then  failure  occurs.  Therefore,  it  is 
important  to  track  the  damage  in  a  structure.  In  the  present 
investigation  some  basic  models  are  proposed  for  the  study  of 
damageable  structure  response.  The  models  are:^Tl)'  a  higher 
order  linear  differential  equation  with  constant  coefficients, 
and  (2)  a  second  order  linear  differential  equation  with  time 
varying  coefficients.  Using  a  digital  computer  a  blast  is  simu¬ 
lated,  and  the  response  of  an  Inelastic  structure  is  computed. 
Noise  signals  are  added  to  these  and  the  results  are  used  to 
simulate  measured  input  and  response.  Next,  using  the  simulated 
input  and  response,  the  parameters  of  the  linear  models  are  iden¬ 
tified  and  the  linear  structure  responses  are  computed.  Mea¬ 
sures  of  these  responses,  Incuding  peak  displacement  and  energy 
dissipated  are  compared  to  the  simulated  response.  It  is  shown 
that  the  models  accurately  simulate  Inelastic  structure 
response.  Moreover,  the  results  of  some  experiments  are 
Included.  The  experiments  show  that  the  energy  dissipated  in  a 
material  specimen  is  related  to  the  damage  level. 
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CHAPTER  1 


1.0  Introduction 

When  a  structure  is  excited  by  an  external  force,  it  exe¬ 
cutes  a  response  determined  by  the  characteristics  of  both  the 
input  and  the  structure.  We  could  predict  the  exact  response  of 
a  structure,  characterized  by  its  geometry  and  its  mutual  proper¬ 
ties,  if  we  could  predict  inputs  exactly;  if  we  had  a  perfect 
model  for  the  structure;  and  if  our  mathematical  computations 
were  correct.  However,  since  inputs  are  random,  we  cannot  per¬ 
fectly  characterize  complex  structures,  and  since  mathematical 
models  are  not  perfect,  we  can  only  estimate  the  response  of  a 
structure. 

In  structural  analysis  we  wish  to  assess  the  response  of  a 
structure  to  dynamic  loads,  such  as  blasts  and  earthquakes.  This 
procedure,  of  course,  requires  the  use  of  a  dynamic  model  which 
will  permit  us  to  predict  the  response  of  a  structure  accu¬ 
rately.  These  structural  models  are  generally  chosen  to  fit 
experimental  data  and  to  simplify  mathematical  computations. 

Most  existing  structures  were  designed  based  on  a  static 
model,  and  although  dynamic  properties  may  be  considered  in  their 
design,  the  designed  parameters  may  be  inadequate  to  predict  the 
response  to  dynamic  load  correctly.  Considerable  work  has  been 
performed  on  identifying  the  parameters  of  mathematical  models 
from  dynamic  experimental  data,  and  various  approaches  have  been 
proposed  for  predicting  system  parameters  based  on  experimental 
data. 

These  identified  parameters  can  be  used  to  predict  the 
dynamic  response  of  a  structure  to  a  different  excitation  than 
that  used  to  test  it.  The  identified  parameters  also  can  be  used 
to  calculate  the  energy  dissipation  in  a  hysteretic  structure 
caused  by  strong  excitation. 


The  energy  dissipated  by  a  structure  during  a  strong  motion 
response  is  an  indication  of  the  structural  damage.  It  is  impor¬ 
tant  to  predict  how  much  damage  occurs  in  a  structure  due  to 
strong  motion  because  the  level  of  damage  is  related  to  the  like¬ 
lihood  of  structural  failure.  When  damage  does  occur,  it  can 
appear  in  different  forms,  such  as  cracks,  permanent  deformation, 
or  change  in  characteristic  frequency. 

While  damage  may  be  located  through  visual  inspection,  such 
an  inspection  may  not  accurately  estimate  either  its  magnitude  or 
exact  location.  Consequently,  it  is  difficult  to  assess  whether 
the  structure  would  survive  normal  design  conditions  or  another 
severe  excitation.  By  visual  inspection,  it  may  be  possible  to 
locate  the  damage  area.  When  damage  occurs  it  is  desirable  to 
estimate  its  magnitude  and  locate  it,  if  possible. 

1.1  Literature  Review 

Part  of  the  energy  dissipated  by  a  structure  is  dissipated 
due  to  hysteretic  behavior  of  the  structural  material.  The  equa¬ 
tion  governing  the  hysteretic  response  of  a  lumped  mass  system  is 
a  second-order,  nonlinear,  ordinary  differential  equation  with 
history-dependent  stiffness  term.  Two  models  which  may  approxi¬ 
mate  the  nonlinear  system  will  be  proposed  in  this  report.  These 
are: 

1.  High-order  equivalent  linear  system; 

2.  Time-varying  parameter  linear  system. 

The  first  model  considered  in  this  paper  is  a  high-order 
equivalent  linear  system.  It  is  assumed  that  the  nonlinear  hys¬ 
teretic  system  is  approximately  governed  by  a  high-order  equiva¬ 
lent  system.  This  model  is  motivated  by  studies  summarized  in 
the  literature.  For  example,  Lutes  and  Hseih  [1]  used  a  third- 
order  linear  system  to  approximate  a  single-degree-of-freedom 
( SOF)  oscillator  with  bilinear  hysteretic  yielding  behavior, 
excited  by  stationary  white  noise.  In  the  linear  system,  certain 


parameters  were  chosen  so  that  the  rms  displacement  and  velocity 
matched  empirical  values  for  the  nonlinear  system.  They  showed 
that  the  third-order  system  gives  a  better  overall  prediction  of 
response  buildup  than  does  either  the  linear  SDF  system  or  a 
two-mode  system. 

Lutes  [2]  used  a  different  type  of  equivalent  linear  system 
to  approximate  the  nonlinear  system.  All  the  methods  Lutes  con¬ 
sidered  defined  the  equivalence  either  in  terms  of  response  dis¬ 
placement  level,  velocity  level,  frequency,  or  a  combination  of 
these.  He  found  that  a  particular  equivalent  linear  system  can 
generally  only  be  expected  to  match  a  limited  number  of  response 
statistics  of  a  particular  nonlinear  system  with  a  particular 
type  of  excitation. 

Wen  [3],  and  Wen  and  Baber  [5,  6]  have  used  the  equivalent 
linearization  method  to  approximately  represent  the  response  of  a 
hysteretic  SDF  system.  They  showed  that  the  third  order,  linear, 
differential  equation  provided  a  satisfactory  representation  of 
the  inelastic,  hysteretic  systems.  This  closed  form  lineariza¬ 
tion  is  relatively  simple  to  formulate  which  allows  ready  exten¬ 
sion  to  multi-degree-of-freedom  (MDF)  systems.  They  showed  that 
the  equivalent  linearization  method  gives  satisfactory  results  at 
all  response  levels  for  response  analysis  of  MDF  deteriorating  or 
non-deteriorating  systems  under  random  excitation. 

Another  study  by  Wafa  [7]  demonstrated  that  the  peak 
response  for  an  hysteretic  SDF  system  excited  by  random  inputs  is 
closely  predicted  by  a  third  order,  linear,  equivalent  system. 
Recent  work  [8]  has  also  shown  that  the  high-order  linear  equiva¬ 
lent  model  provides  a  good  approximation  to  the  hysteretic  system 
when  the  energy  dissipated  and  frequency  shift  are  concerns. 
Significantly,  the  results  established  that  the  parameters  of  a 
higher  order  system  can  be  identified  by  using  a  frequency  domain 
method  even  when  noise  is  present  both  in  the  forcing  and 
response  signals.  In  contrast,  the  time  domain  approach  yields 
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poor  results  in  the  presence  of  noise.  Because  of  the  frequency 
domain's  preferable  application,  it  will  be  used  to  do  the 
analysis. 

The  second  model  is  motivated  by  the  fact  that  structures 
may  exhibit  time-variant  nonlinear  response  to  strong  motion. 

This  implies  that  structure  deterioration  was  in  progress  during 
the  large  amplitude  motion.  Such  a  phenomenon  has  been  recog¬ 
nized  and  studied  in  the  past.  For  example,  Udwadia  and  Trifunac 
[9]  and  Iemura  and  Jennings  [10]  carried  out  an  analysis  to  char¬ 
acterize  such  behavior  in  terms  of  a  quasi -time-variant  linear 
formulation  based  on  the  data  obtained  from  the  San  Fernando 
earthquake  of  February  9,  1971. 

In  another  study  Townsend  and  Hanson  [11]  demonstrate  time- 
varying  hysteretic  loops  by  the  experimental  test  of  reinforced 
concrete  beam-column  and  T-shaped  specimens  under  different  load¬ 
ing  conditions.  In  addition,  Uzumeri  [12]  has  also  shown  the 
same  behavior  for  an  experimental  study  of  cast-in-place  rein¬ 
forced  concrete  beam-column  joints  subjected  to  simulated  seismic 
loading. 

Based  on  the  above  referenced  investigations  involving 
time-varying  parameters  system,  we  anticipate  that  the  time-vary¬ 
ing  parameter  model  will  provide  a  good  representation  of  a  hys¬ 
teretic  system. 

1.2  Objective 

The  determination  of  the  system  parameters  from  suitable 
experimental  observations  is  a  fundamental  problem  in  engineer¬ 
ing.  Obtaining  a  good  representation  of  a  system  requires  all 
the  proper  information,  such  as  well  measured  data  and  a  suitable 
model. 


The  objective  of  this  study  is  to  justify  two  possible  mod¬ 
els  to  characterize  the  behavior  of  a  system.  The  relative  mer¬ 
its  of  each  model  are  discussed.  Extensive  numerical  experimen¬ 
tation  using  simulated  data  is  also  presented  in  order  to  inves¬ 
tigate  their  feasibility  and  accuracy.  This  study  will  demon¬ 
strate  how  well  the  system  parameters  can  be  identified  with  and 
without  noise  in  the  measurements.  Once  the  parameters  are 
known,  the  energy  dissipated  by  the  system  can  be  computed. 

Based  on  the  computed  results,  one  can  compare  how  well  the  mod¬ 
els  performed  for  a  given  set  of  data.  The  ultimate  goal  of  this 
study  is  to  establish  structural  models  useful  for  other  pur¬ 
poses,  such  as  prediction,  design,  control,  and  damage 
assessment. 

The  present  research  has  been  aimed  at  the  analysis  of  dam¬ 
age  accumulation  in  concrete  structures.  It  is  assumed  that  when 
a  concrete  structure  dissipates  energy,  it  accumulates  damage. 

To  justify  this  assumption  some  physical  experiments  have  been 
performed.  Specifically,  concrete  cylinders  have  been  subjected 
to  cyclic  loading.  The  energy  dissipated  in  each  cylinder  was 
measured  and  the  level  of  residual  strength  in  each  cylinder  was 
determined  after  the  load  cycling  was  completed.  The  residual 
strength  was  plotted  versus  energy  dissipated.  When  the  reduc¬ 
tion  in  strength  is  taken  as  a  measure  of  damage,  this  plot 
reveals  the  damage  caused  by  energy  dissipation. 


CHAPTER  2 

HIGH  ORDER  EQUIVALENT  LINEARIZATION 

2.1  Model 

The  differential  equation  governing  the  response  of  a 
single-degree-of-freedom  (SOF)  system  is 

mz  +  u  =  f  (2-1) 

where  m  is  the  mass  of  the  structure,  f  is  the  forcing  function, 
z  is  the  displacement  response,  dots  denote  differentiation  with 
respect  to  time,  and  u  is  the  restoring  force  of  the  structure. 
Equation  (2-1)  can  be  used  to  model  the  actual  system  in  which  u 
can  be  a  very  complicated  function.  In  the  present  study,  the 
hysteretic  restoring  force,  u,  is  modeled  by  using  the  equation 

fa  v<J)  =  w  * 2  (2'2> 

where  the  Cj,  j  =0,  1,  ...,  M+l  are  the  constants  governing 
the  system  restoring  force  characteristics,  where  u(J)  denotes 
the  jth  time  derivative  of  u,  and  M  is  a  constant  denoting  the 
order  of  approximation  provided  by  the  linear  system.  The  reason 
for  using  this  model  to  represent  the  hysteretics  system  is  that 
it  displays  a  hysteretic  character  that  can  be  made  to  match  the 
character  of  an  inelastic  structure. 

Consider  the  case  where  M  is  equal  to  0.  The  model  in  Equa¬ 
tions  2.1  and  2.2  becomes 

mz  +■£■!■  z  +  z  s  f  (2-3) 

co  co 
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DISPLACEMENT 

Figure  2.1  Total  restoring  force  versus  displacement  for  second-order  system 


This  is  simply  the  second-order  linear  differential  equation 
governing  the  SDF  system.  However,  when  the  response  of  the 
actual  system  is  linear  and  damping  is  viscous,  the  model  of 
Equation  (2-3)  represents  the  actual  system.  The  restoring  force 
function  for  this  system  is  u  =  (ci/co)z  +  (l/co)z.  This  model 
displays  the  hysteretic  behavior  as  shown  in  Figure  (2.1). 

When  the  constant  M  is  chosen  as  1  in  equation  2-2,  the 
model  becomes 


mz  +  u  =  f 

I  I 

CoU  +  ClU  =  C2Z  +  Z 


(2-4) 


Combining  these  equations  results  in 


T-i-  mz  +  mz+^z+-iz  =  f+  4if 


(2-5) 


The  parameters  of  the  system  in  Equation  (2-5)  can  be  chosen  so 
that  the  model  represents  the  hysteretic  system  as  well  as 
possible. 

For  example.  Figure  (2.2)  shows  the  hysteretic  properties 
for  an  SDF  system  by  plotting  the  restoring  force  versus  dis¬ 
placement.  The  parameters  for  the  system  and  the  forcing  input 
are  given  in  Figures  (2.2).  This  study  will  consider  both  the 
equations  (2-3)  and  (2-5). 

The  parameters  can  be  identified  by  using  the  least-square 
identification  criterion.  Since  most  observed  data  include  a 
certain  percentage  of  noise,  the  frequency  domain  approach  will 
be  used  to  perform  the  parameter  identification.  However,  a 
practical  problem  of  accuracy  arises  in  the  parameter  identifica¬ 
tion  problem  for  high-order  linear  systems.  Therefore,  the 


MV... 
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high-order  systems,  for  M  >  1,  will  not  be  considered  in  the 
following  discussion.  It  is  anticipated  that  in  some  cases,  the 
third-order  linear  approximation  (Equation  2-5)  provides  a 
better  representation  of  the  hysteretic  system  than  the  second- 
order  linear  system  (Equation  2-3). 

2.2  Identification  of  Parameters 

The  problem  of  parameter  identification  can  be  posed  as  one 
class  in  the  broader  topic  of  optimization.  The  object  of  param¬ 
eter  identification  is  to  make  inferences  about  the  real  world 
and  mathematical  models  on  the  basis  of  measured  input  data.  The 
measured  data  in  this  study  were  assumed  available,  and  were 
simulated  to  represent  the  field  data. 

First,  it  is  assumed  that  there  is  no  noise  present  in  the 
measured  data.  Then,  noise  data  are  introduced.  Note  that  the 
measured  response  data  are  given  as  acceleration  values.  This  is 
realistic  since  the  structural  response  acceleration  is  often  the 
measured  quantity  in  an  experimental  test. 

2.2a  Second  Order  System 

Now,  consider  the  second-order  model.  Equation  2-3,  Equa¬ 
tion  2-3  can  be  simplified  by  taking 

<£■*•■  •>■!}  <2-6> 

The  equation  governing  motion  of  the  system  becomes 

if  I 

mz  +  aiz  +  aoz  =  f  (2-7) 

Fourier  transform  both  sides  to  obtain 

m(iw)2Z(u>)  +  ai(iu)Z(w)  +  aoZ(u)  =  F(o>)  (2-8) 

where 


in 


Z(u)  =  /  z(t)  e"1u>t  dt  -»  <  a»  <  » 
.« 

ao 

F(u)  =  /  f(t)  e"la,t  dt  -»  <  u  <  « 
-00 


(2-9) 


are  the  Fourier  transforms  of  z(t)  and  f(t),  respectively.  This 
equation  can  be  rearranged  and  combined  with  Z(u)  and  F(u)  terms 
on  one  side  of  the  equation  to  obtain 

2  F 

(-mu  +  ao  +  aiiu)  =  ^  (2-10) 


Multiply  each  side  of  the  equation  by  its  complex  conjugate  to 
obtain  the  modulus  squared. 


l-m<u^  +  ao  +  ai  iul 

1  1  izr 


(2-11) 


Evaluate  the  left  side  and  let  |F|^/|Z|2  equal  Q(u)  to  obtain 

(2-12) 


o  u  o 

(ao  -  mu  )  +  (aiu)4-  =  Q(u) 


This  equation  can  be  used  to  identify  the  parameters  of  second- 
order  linear  system.  However,  this  equation  is  exactly  satisfied 
if  and  only  if:  1)  the  system  under  consideration  is  linear;  2) 
all  measurements  are  noise  free;  and,  3)  the  Fourier  transforms 

II 

of  z(t)  and  f(t)  used  In  Equation  (2-12)  are  exact.  When  these 
requirements  are  not  satisfied,  Equation  2-12  will  include  an 
error  term  (or  noise  term).  This  practical  case  is  usually  the 
one  that  needs  the  most  attention. 

When  noise  is  present  on  the  measured  input  and  response. 
Equation  2.12  can  be  written  as 


2  6  2 

(ao  -  mu  )  +  (aiu)  =  Q(u)  +  e(u) 


(2-13) 
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Note  from  Equation  2-7  that  ao  is  the  equivalent  stiffness  and 
aj  is  the  equivalent  damping  for  the  second-order  linear  system. 
Therefore,  ao  is  greater  in  magnitude  than  aj .  In  view  of  this 
and  the  form  of  Equation  2-13,  a0  can  be  estimated  by  noting  the 
frequency  where  Q(u)  +  e(u>)  is  a  minimum  whenever  the  equivalent 
damping  factor  is  much  less  than  1  (say  less  than  0.2).  This 
will  be  true  in  most  civil  engineering  systems. 

Denote  the  frequency  where  Q( oj)  +  e( m)  is  a  minimum  by 
ofo.  Equation  2-13  shows  that,  approximately, 

a0  =  mu)m2  (2-14) 

since  the  first  term  on  the  left  is  approximately  zero  when 
Q(u>)  +  e(oj)  is  a  minimum.  Substitute  Equation  2-14  into  Equation 
2-13;  this  yields 

m2(<i^2  -  u2)  +  (ai<u)2  =  Q(u»)  +  e(oi)  (2-15) 

Now,  it  is  necessary  to  find  the  coefficient  ax  which  minimizes 
the  e(oj). 

The  coefficients  ax  can  be  evaluated  using  a  least-squares 
approach,  where  the  integral  of  e^(u>)  over  a  specific  range  of 
frequencies  is  minimized.  Based  on  Equation  2-15,  set 


where  a>a  amd  qj  are  lower  and  upper  bound  frequencies, 
respectively.  This  frequency  band  is  chosen  so  that  the  system 
response  behavior  can  be  fully  characterized.  It  is  anticipated 


K2- 


„2> 


-  Q(u>)  I  du>  = 


e2(u>) 


(2-16) 
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that  the  frequency  band  includes  the  natural  frequency  for  a 
linear  or  slightly  non-linear  system.  For  a  highly  nonlinear 
system,  the  characteristic  frequency,  u^,  will  shift.  However, 
the  frequency  band  can  be  located  by  finding  the  frequency  where 
Q(w)  +  e(w)  is  a  minimum,  and  selecting  the  frequency  band  around 
this  frequency. 

In  general,  the  lower  frequency,  wa,  is  located  at  a  point 
where  its  corresponding  Q{wa)  +  e(“a)  value  is  about  5  times 
as  great  as  the  minimum  value  of  Q(w)  +  e(u>);  and  is  the 
higher  frequency  where  Q(^b)  +  e(^)  is  about  5  times  greater 
than  the  minimum  value  of  Q(w)  +  e(to).  This  method  is  used  to 
select  u>a  and  based  on  an  approximate  linear  analysis. 

When  u>a  and  are  chosen  in  this  manner,  the  interval  (wa, 

Mb)  will  be  approximately  the  half  power  bandwidth  of  the  sys¬ 
tem.  This  frequency  interval  reflects  the  characteristics  of  the 
system. 

Take  the  first  partial  derivative  of  e2  with  respect  to 
ai2,  and  set  it  equal  to  0.  The  result  is 


&  ■  -n 


ai2u>2  +  m2  (*/  -  u2)  -  Q(< 


du>  =  0  (2-17) 


Simplify  this  to  obtain 


w  “h 

j  »4  d.  -  f  (q 


2  ,  2 
m  ( w 


J)2\  J  dw  (2-17a) 


Integrate  the  equation  where  possible  to  get 


%  =  qb  % 

where  qa  is  a  coefficient  less  than  one  and  is  a 
coefficient  greater  than  one.  This  equation  can  be  further 
simpl if ied 


This  equation  provides  the  value  of  aj.  This  is  the  best  esti¬ 
mator  in  the  least-squares  sense.  All  the  parameters  in  this 
equation  are  known  except  the  integral  of  the  oj^Q(to)  term  which 
can  be  evaluated  numerically. 

2.2b  Third  Order  Equation 

The  second-order  linear  ordinary  differential  equation  may 
not  be  considered  an  accurate  representation  of  the  hysteretic 
system.  It  is  hoped  that  the  third-order  linear  system  may 
improve  the  accuracy  in  some  sense. 

Equation  2-5  can  be  simplified  by  taking 
1  c2  ci 

a°’  c7  '  *2  (M1) 
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Then  Equation  2-5  becomes 

in  H  i  i 

ma2  z  +  mz  +  a^  +  a0z  =  f  +  a2f 
Fourier  transform  both  sides  to  get 

(ma2  (iu)3  +  m(iio)2  +  a,  (iu)  +  a0 


|1  +  a2  (iaj 

The  symbols  used  in  this  equation  have  the  same  meaning  as  in 
earlier  equations.  Multiply  each  side  of  the  equation  by  its 
complex  conjugate  to  yield  the  modulus  squared 


3  2  ^ 

ima2uj  -  mw  +  ia^oi  +  a2|  |p|2 

1  +  ia2  u)|2  |Z|2 


(2-24) 


Evaluate  the  left  hand  side  and  let  |F|2y/jZ|2  be  replaced  by 
0(<u)  to  obtain  ; 


(a0  -  moi2)  +  (»)2  ( a x  -  mco2a2)2 
- - - 

1  +  (a2(u) 


=  Q  (w) 


(2-25) 


This  equation  governs  a  third-order  linear  system  in  the  fre¬ 
quency  domain.  Measured  data  can  satisfy  this  equation  exactly 
if  and  only  if:  1)  the  system  under  consideration  is  linear;  2) 
all  measurements  are  noise  free;  and,  3)  the  Fourier  transforms 
used  to  define  Q(w)  are  exact.  These  conditions,  however,  are 
not  usually  met.  In  fact,  the  purpose  of  this  investigation  is 
to  use  the  higher-order  linear  system  to  represent  an  hysteretic 
system.  Therefore,  measured  data  do  not  usually  satisfy  the 


above  equation.  To  account  for  this  explicitly.  Equation  2-25  is 
written 


2  2  2  p  ^ 

(ao  -  mw  )  +  w' ( a i  -  mw  a2) 

1  +  (a2“)^ 


Q(  o>)  +  e(w) 


(2-26) 


e(w)  is  a  noise  term  which  must  be  minimized  by  the  proper  choice 
of  system  parameters.  This  equation  can  be  used  to  identify  the 
system  parameters  following  several  approaches.  Two  of  these  are 
summarized  below.  One  approach  approximates  certain  terms  in 
Equation  2-25  to  obtain  estimates  for  the  system  parameters, 
while  the  other  approach  uses  a  search  technique  to  estimate 
system  parameters. 

The  first  method  to  be  investigated  is  an  approximate  tech¬ 
nique.  When  a2  is  small  compared  to  the  characteristic  fre¬ 
quency,  wm  of  an  SDF  system,  the  (a2w)^  term  in  the  denomi¬ 
nator  can  be  neglected.  This  is  usually  true  when  non  linear 
deformation  is  not  too  large.  Eliminate  the  (a2w)2  term  in 
Equation  2-26  to  obtain 


2  2 

(ao  -  mw2)  +  w2  (ai  -  mw2a2)  =  Q(w)  +  e(w)  (2-27) 

When  ai  and  a2  are  small  compared  to  ao  (which  is  usually  the 
case,  the  minimum  of  the  left-hand  side  occurs  near  the  frequency 


w  = 


(2-28) 


As  previously  described,  the  characteristic  frequency,  w^,  can 
be  found  vtfien  the  Q(w)  +  e(w)  term  is  a  minimum.  In  terms  of 
wm,  ao  can  be  written 

ao  =  mw^2  (2-29) 
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Expand  the  second  term  on  the  left  side  of  Equation  2-27  and  use 

2  4  2 

the  result  of  Equation  2-29.  Neglect  the  nun  dj  term.  Then 
Equation  2-27  becomes 


m2  (iom2  -  u^2)  +  w2  {a*2  -  2ma)2a|a2)  =  Q(oj)  +  e(w)  (2-30) 


Let  a^  =  bi  and  aia2  =  b2,  then 


(2-31) 


oj2  (bi  -  Sn to2 b 2 )  Q(oi)  -  m2(<^2  -  w2)  j=  e(w)  (2-32) 

When  this  expression  is  evaluated  at  the  discrete  frequency 


oj  =  a*,  the  result  is 


\  <bl  -  2nl“k2b2>  -(^k  '  m2(“m2  '  “k2)  )  -  El 


(2-33) 


The  system  parameters,  aj  and  a2,  can  now  be  identified  as 
those  which  minimize  the  sum  of  the  squares  of  the  terms  in 
Equation  2-33.  Consider  a  sequence  of  discrete  frequencies 
uniformly  spaced  in  the  interval  (wa,  u^).  These  are  <ua, 
ua  +  Aid,  o)a  +  2auj,  etc.  Define 


{b}  =  (bx  b2)T 


(2-34a) 


-  2m a). 


w_  +  Ab>  -2m(w.  +  Au) 

a  a 


fx,l  =  <*»  +  2Aai 


2m(w,  +  2au>) 


(2-34b) 


(2-34c) 


by 


Aa 


-  m 


!(  (»m2-»a2)f 


{*?}  = 


'(ua+A<j 


-  m  \  m 


(“m  •  (uja+&uj)< 


)' 


^o»  +  2Aw 

a 


m 


(“m2  ‘  <“a  *  2s“)2) 


to. 


2  /  2  2v ‘ 

-  m  -  tob  ) 


M  =  (£(U  c_ 

a  a 


c*  i  •  •  •  *  •  6  J 

u,  u>,+Auj  io^' 


(2-34d) 


Note  Aw  is  equal  to  ( 2-rr/T)  and  T  is  the  total  duration  of  the 
excitation.  (aia,  defines  the  range  of  frequencies  over  which 
the  system  is  analyzed.  As  in  the  identification  of  parameters 
of  the  second-order  linear  system,  only  a  portion  of  the 
frequency  range  is  used  in  the  parameter  identification. 

This  frequency  range  can  be  chosen  as  before.  In  terms  of 
the  matrices  defined  above.  Equation  2-32  can  be  written  at 
discrete  frequencies  as 


Cxx3  lbl  -  {x2 1  =  M 


(2-35) 


The  vector  {b}  can  be  found  by  minimizing  the  sum  of  the  squares 
of  {e}.  This  is  =  {e}^  {e}.  The  vector  {b}  which 


minimizes  is 


fb}  ■  <tx13Ttx13)‘1Cx1]T  (x2) 


(2-36) 


When  bi  and  b2  have  been  computed,  ai  and  a2  can  be  found  from 
Equation  2-31. 


a i  =  /57,  &2  =  b2/ a i 


(2-37) 
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This  solution  provides  approximate  values  for  the  parameters  in 
the  third-order  linear  system  which  represents  the  hysteretic 
system. 

Another  approach  can  be  used  for  the  estimation  of  param¬ 
eters  in  the  higher  order  linear  system.  This  is  a  search  pro¬ 
cedure  which  iteratively  estimates  the  parameter  values.  Equa¬ 
tion  2-26  can  be  rewritten 


2  2 
(a0  -  m<u  )  +oj  (a^  -  mw  32) 

«(«)  »  Q(»)  -  - 7- -  (2-38) 

1  +  ( a£U) ) 

This  quantity  is  a  measure  of  the  mismatch  between  the  measured 
data,  reflected  in  Q(w),  and  the  model,  reflected  in  the  second 
term  on  the  right  side  of  Equation  2-38. 

This  mismatch  can  be  either  positive  or  negative  and  can  be 
used  to  define  one  measure  of  the  difference  between  the  model 
and  the  measured  data  over  a  range  of  frequencies.  This  measure 
is 


2 


e 


e(«Jk) 


(2-39) 


where  the  sun  is  taken  over  those  discrete  frequencies  in  the 
interval  (wa,  0^).  This  is  the  square  error  of  the  model. 

This  error  is  minimized,  however,  when  the  model  parameters  are 
chosen  to  satisfy  the  sequence  of  equations 


2  2  2 
3e_  Q  _  3e*  =  3e 

3a0  3ai  882 


(2-40) 


The  parameters  ao,  ai,  and  a2,  which  satisfy  these  equations 
establish  a  model  which  is  optimal  in  a  least  squares  sense. 


Equations  2-40  can  be  solved  numerically  using  a  search 
technique.  A  computer  program  was  written  to  solve  Equation 


2-40.  The  program,  included  in  the  appendix,  uses  Newton's  meth 
od  to  search  for  the  solution.  The  analysis  procedure  followed 
in  the  computer  program  is  identical  to  that  used  in  solution  of 
the  problem  summarized  in  the  following  section.  The  steps  in 
the  solution  procedure  are  listed  at  the  end  of  Chapter  3. 


CHAPTER  3 


SYSTEM  WITH  TIME-VARYING  PARAMETERS 

3.0  Time-Varying  Parameters  Model 

Sructures  may  exhibit  time-variant  nonlinear  response  to 
strong  motion  excitation.  Time-varying  structural  properties 
were  not  considered  in  the  previous  chapters.  In  this  chapter,  a 
structure  is  modeled  as  a  time  variant  si ngle-degree-of -freedom 
(SOF)  oscillator,  and  a  methodology  is  introduced  to  determine 
its  parameters  using  the  observed  data.  It  is  important  to 
introduce  a  technique  which  can  be  applied  when  noise  is  present 
in  the  measured  data. 

To  demonstrate  this  procedure,  consider  an  SDF  linear  system 
with  mass  m.  Let  the  damping  and  stiffness  parameters  for  this 
system  be  time  varying.  Its  equation  of  motion  is 

mz  +  C(t)z  +  K(t)  z  *  f  (3-1) 

in  which  z  is  the  displacement  response  of  the  system;  C(t)  and 
K(t)  are  time  variant  damping  and  stiffness  function  of  the  sys¬ 
tem,  respectively;  and  f  is  the  forcing  funtion.  It  is  proposed 
that  this  equation  be  used  to  model  the  behavior  of  a  system 
governed  by  Equation  2-1.  Observe  the  system  from  time  0  to  T 

I 

and  assume  that  z(0)  =  0  and  z(0)  =  0. 

The  functions  (C(t)  and  K(t)  are  assumed  to  have  the  form 

C(t)  =  (1  +  at)  co 

K(t)  =  (1  +  3t)  ko  (3-2) 


Here,  Co  and  ko  are  the  damping  constant  and  the  stiffness  con¬ 
stant,  respectively,  a  and  6  are  consant  coefficients  which  are 
usually  much  less  than  one. 


In  many  practical  cases,  it  is  observed  that  the  structure 
displays  an  increase  in  damping  and  a  decrease  in  stiffness  when 
the  structure  excites  an  inelastic  response.  This  implies  a  is  a 
positive  constant  and  S  is  a  negative  constant. 

In  this  study,  although  a  and  B  will  be  considered  as  small 
values,  they  will  be  large  enough  to  influence  the  system's  prop¬ 
erties.  This  permits  treatment  of  Equation  3-1  as  a  perturbed 
differential  equation.  When  a  and  B  are  both  equal  to  0  (unper¬ 
turbed),  the  equation  3-1  is  simply  a  second-order  differential 
equation  with  constant  coefficients  which  can  be  easily  solved. 

The  solution  of  Equation  3-1  can  be  written  in  the  form  (for 
example  reference  18) 

z  ■  Zo  +  <*  za  +  B  Zg  +  high-order  terms  (3-3) 

Since  a  and  B  are  small,  the  high-order  terms  will  be  neglected. 
Substituting  Equation  3-3  into  3-1  and  expanding  yields 

m(zo  +  azQ  +  3z0)  +  co  (1  +  at)  (zo  +  azQ  +  BZg) 

+  ko  (1  +  3t)  (zo  +  aza  +  Bz0)  =  f  (3-4) 

Moving  the  force  term  to  the  left  side  of  the  equation,  grouping 
coefficients  of  the  terms,  1,  a,  and  B,  then  equating  the  coeffi¬ 
cients  to  zero  results  in 

ll  I  . 

mzo  +  cozo  +  kozo  =  f  ( 3— 5a) 

Ml  • 

mza  +  c0za  +  koza  =  -  cotzo  (3-5b) 

Ml 

mzB  +  cozB  +  koZg  =  -  kotzo  (3-5c) 

These  equations  approximately  govern  the  system's  response  when 
time  variation  of  the  parameters  is  linear,  as  shown  in  Equation 
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3-2.  If  the  excitation  and  the  system's  response  are  known,  then 
these  equations  can  be  used  with  a  time  domain  parameter 
identification  procedure  to  estimate  the  system  parameters. 
However,  when  a  time-domain  parameter  identification  approach  is 
used,  problems  arise  if  noise  is  present  in  the  measured  input 
and  response  signals  (see  Reference  8). 

A  frequency  domain  approach  to  the  identification  of  system 

i 

parameters  is  pursued.  Therefore,  the  equations  of  motion  are 
transformed  to  the  frequency  domain. 


qo  # 

Z0(oj)  =J  z 0(t )  e~mt  dt 


(3-6a) 


Z0M  =/  2Q(t)  e  dt 


Z0U)  =  f  Zg(t )  e"1a)t  dt 

MOO 

define  the  Fou/ier  transforms  of  z0(t),  za(t),  and  Zg(t). 


(3-6b) 


(3-6c) 


And  let 


z°(t)  s  "j—  f  Z0(.)  ela)t  dm 

SmOO 

Z*M  =  27  /  Za^u)  eiWt  dw 

MOO 


( 3-7  a) 


(3-7b) 


(3-7c) 
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define  the  inverse  Fourier  transforms.  Then  the  Fourier  trans¬ 
form  of  Equation  3-3  is 

Z(“)  =  Zo  +  aZa  +  ez3  (3-8) 

It  can  be  shown  that  the  Fourier  transform  of  Equations  3-5a 
through  3-5c  are  given  by 

-nto^Zo(w)  +  co  iwZo(w)  +  ko  Z o ( a>)  =  F(w)  (3-9a) 

-m“^Za  ( “>)  +  co  i“Za(<»>)  +  koZa(u)  =  co^Zo(w)  +  <*>Z'o(“)) 

(3-9b) 

O 

-m<*>  Z3(“)  +  co  i“Z3(w)  +  koZ3(w)  =  -i  ko  Z'o(<*>)  (3-9c) 


Now  solve  the  Equations  3-9  simultaneously.  The  result  is 

Zo(«)  =  H(o»)  F(oi)  (3-10a) 

z0(»)  *  CoH(oj)  |h(u)F(u»)  +  u(h'(w)F(uj)  +  H( u>) F ' ( u>)  (3-10b) 


Z3(«)  =  -ikoH(w)  (H*F(u>)  +  H(»)  F*(u»)) 


(3- 10c) 


where  H(a>)  is  frequency  response  function  and  H * ( u>)  is  its  first 
derivative.  These  can  be  written  in  the  forms 


H(o>)  =  [(ko  -  mor)  +  i(wco)] 


-1 


H'(w)  *  (2ma>  -  ic o)  C ( ko  -  mar)  +  i(wco)] 


F(u>)  is  the  Fourier  transform  of  f(t)  and  F'(u)  is  its 
derivative. 


oo 

F(w)  -J  f(t)  e 


-iwt 


dt 


(3-lla) 
( 3- lib ) 


(3-12a) 


OA 


dt 


( 3-12b) 


oo 

F1  (to)  =  -i  J  tf(t)  e"lu)t 

a  00 

Substitute  the  results  from  Equations  3-10  into  Equation 
3-8.  This  yields 

Z( oj)  =  H(uj)F(oj)  +  aCo  H(oj)  |h(oj)F(u)  +  <d(H' (u)F(a>) 

+  H(u)F'(u))  +  ( - i k O0 )  H{u)  (  H*U)F(u>)  +  H(u)F'(u))  (3-13) 

This  is  the  approximate  frequency  domain  expression  for  the  solu¬ 
tion  of  Equation  3-1.  It  is  considered  accurate  when  both  a  and 
g  are  small.  The  displacement  response  also  can  be  obtained  by 
inverse  Fourier  transformation  of  the  Equation  3-13.  Equation 
3-13  is  used  in  the  identification  process.  Its  use  will  finally 
lead  to  the  estimation  of  the  parameters  from  a  sequence  of  mea¬ 
sured  data. 

3.1  Identification  Procedure 

The  method  described  above  provides  the  solution  for  the 
Equation  3-1  in  the  frequency  domain.  When  the  measured  values 


of  f(t)  are  used  to  estimate  FU)  and  the  result  is  used  in  Equa¬ 
tion  3-13  to  obtain  Z(oj),  this  Z(a>)  will  not,  in  general,  match 
the  Z( to)  estimated  from  the  measured  Z(t).  Moreover,  the 
calculated  |Z(w)|  will  not  match  the | Z(oj)  |  obtained  from  measure¬ 
ment.  The  reasons  for  this  mismatch  are  that  (1)  noise  is  inevi¬ 
tably  present  in  the  measured  input  and  response.  (2)  the  mathe 
matical  model  is  linear,  yet  the  measured  data  come  from  non¬ 
linear  structures  and  (3)  the  discrete  Fourier  transform  of  a 
time  series  is  used  to  represent  the  continuous  Fourier  trans¬ 
form.  In  the  following,  a  brief  theoretical  background  is  pre¬ 
sented  together  with  the  simple  description  of  the  procedure  for 
finding  the  unknown  parameters. 


An  equation  defining  the  mismatch  between  the  measured  data 
and  the  model  of  Equation  3-1  can  be  established.  Let 
|z(m)(a))|  be  the  modulus  of  the  Fourier  transform  of  the 
measured  structural  response  data.  Let  |Z(w)|  be  the  modulus  of 
the  function  obtained  when  the  Fourier  transform  of  the  measured 
input  data  is  used  in  Equation  3-13.  The  difference  between 
these  function  is  defined 

e(u>)  =  |Z(  w)  j  -  |Z(u»)(m)|  (3-14) 

When  the  discrete  Fourier  transform  is  used  to  approximate 
the  continuous  Fourier  transform  of  a  measured  or  theoretical 
signal,  it  is  defined  at  a  discrete  set  of  frequencies, 
a*  =  kAw,  k  =  0,1,..., n-1.  Here  n  and  Aw  relate  to  the  time 
signal  z(t)  and  its  discretization.  It  is  assumed  that  z(t)  is 
available  on  the  interval  (0,  T)  and  is  represented  by  the 
discrete  set  of  values  Zj,  j=0,.../n-l.  Thus,  n  is  the  number 
of  points  where  the  signal  is  represented.  Aw  is  given  by  2ir/T. 
At  a  particular  frequency  w  =  wk,  Equation  3-14  becomes 

ek=  |Zk(w)|  -  |Zk(m)(w)|  (3-15) 

ek  can  be  positive  or  negative.  A  quantity  which  is 
always  non-negative  and  which  summarizes  the  differences  between 
the  measured,  )Zkm^  |,  and  the  theoretical,  |Zk|,  structural 
responses  in  the  frequency  domain  over  a  range  of  frequencies  is 
given  by 


This  is  the  square  error  between  measured  data  and  the 
model.  This  error  can  be  minimized  by  properly  choosing  the 
parameters,  k0,  Co,  a,  and  6.  A  method  that  chooses  the 
parameters  this  way  is  a  least  squares  method. 


The  range  of  index  values,  k,  over  which  the  above  sum  is 
taken,  is  not  specified  in  Equation  3-16.  Equation  3-16  need  not 
be  summed  from  0  to  n.  Rather,  the  summation  should  be  carried 
out  over  the  range  of  freuencies  which  includes  those  values  of 
Z|<  containing  significant  information  on  the  behavior  of  the 
system.  In  general,  this  is  the  band  of  frequencies  surrounding 
the  characteristic  frequency  of  the  system. 

Now,  one  can  choose  k0,  c0,  a,  and  0,  as  those  constants 
which  satisfy 


2  2  2  2 
3e  =  i£_  =  l£_  =  l£_  = 

3k  o  3C  o  3a  30 


(3-17) 


The  ko,  Co,  a,  and  0,  can  be  located  using  a  search  technique. 

To  simplify  the  analysis,  Newton's  method  is  used  to  minimize 
e2  with  respect  to  ko,  Co,  a,  and  0. 

Newton's  method  converges  very  rapidly  once  an  iterate  is 
fairly  close  to  the  solution.  The  formal  simplicity  and  its 
great  speed  are  the  reasons  why  Newton's  method  is  used  in  this 
study. 

To  assure  convergence  in  the  numerical  analysis,  it  is 
important  to  choose  the  initial  iterate  properly.  A  more 
detailed  discussion  of  the  numerical  procedures  will  be  presented 
later,  in  the  numerical  examples. 

The  steps  in  the  numerical  analysis  are  as  follows: 

1.  Make  the  initial  guesses  at  the  parameter  values,  k0, 
c  o ,  a,  and  0. 

2.  Choose  the  computation  increments  Ako,  acq,  Aa,  and  A0. 


3.  Choose  the  desired  accuracy  measure  (used  to  judge 
convergence) . 


4.  Compute  the  partial  first  and  second  derivatives  of 
with  respect  to  k0  using  central  difference  formulas. 

5.  Use  Newton's  method  to  minimize  with  respect  to 
k0. 

6.  Repeat  steps  4  and  5,  this  time  minimizing  with  respect 
to  Co,  then  o,  then  $. 

7.  Check  the  result  to  convergence. 

a.  If  convergence  has  occurred,  then  stop  the 
analysis. 

b.  If  convergence  has  not  occurred,  then  repeat  steps  4 
through  6. 

A  computer  program  to  execute  the  procedure  described  above 
has  been  written.  This  proqram  is  named  PUR  and  a  listing  is 
included  in  the  Appendix. 


CHAPTER  4 


EXAMPLES 

4.0  Numerical  Examples 

In  this  chapter,  numerical  examples  are  presented  to  demon¬ 
strate  the  use  of  the  previously  described  analytical  proce¬ 
dures.  Two  examples  demonstrating  the  frequency  domain  approach 
to  parameter  identification  are  summarized.  Both  solutions 
employ  the  sane  excitation  input.  The  response  signals  identi¬ 
fied  were  drawn  from  two  sources.  One  source  was  a  bilinear  hys- 
teretic  response  while  the  other  was  a  time  varying  linear  sys¬ 
tem.  Use  of  these  two  sources  may  provide  improved  understanding 
of  the  feasibility  of  the  identification  procedure  for  specific 
measurements.  The  examples  show  the  identifications  of  the 
parameters  for  linear  and  hysteretic,  single  degree  of  freedom 
structures,  when  measurement  noise  is  and  is  not  present. 

The  input  used  to  excite  the  SDF  system  in  all  numerical 
examples  is  a  decaying  exponential,  oscillatory  function.  It  is 
generated  using  the  formula 


Where  a,  c.,  j=l,...,N,  and  <u.,  j=l,...,N  are  constants. 

J  J  J 

j*l,...N,  are  phase  angles  which  are  random  variable  realiza¬ 
tions;  these  random  variables  are  independent  and  uniformly 
distributed  on  the  interval  (0,  2tt).  a  is  a  decay  rate.  The 
Cj,  j=0,...N,  are  constants  which  determine  the  amplitudes  of 
the  excitation.  All  the  values  of  Cj  are  taken  as  equal  to  c 
in  all  cases  for  the  examples,  ujj,  j=l,...N,  are  the  frequen¬ 
cies  where  the  excitation  has  power,  uj,  j=l,...,N,  are 
equally  spaced  in  the  interval  including  the  characteristic 
frequency  of  the  system  being  analyzed. 

The  forcing  function  defined  above  was  generated  at  discrete 
times.  Specifically,  f(t)  was  evaluated  at  the  times  t=t.=tAt, 


1=0,..., N-l.  A  computer  program,  named  FORCE,  which  generates 
the  excitation  of  Equation  4-1  was  used  in  these  numerical 
examples. 

Two  distinct  signal  types  were  identified  in  the  numerical 
examples.  The  first  example  used  a  computer  program,  named 
BILIN,  to  compute  linear  and  nonlinear  response.  BILIN  can  be 
used  to  find  the  displacement,  velocity,  and  acceleration 
response  of  a  given  bilinear  hysteretic  system  to  an  arbitrary 
input.  It  also  computes  the  energy  dissipated  by  the  structure 
during  the  response.  The  second  example  used  a  computer  program, 
named  TIMEVA,  to  compute  the  response.  This  program  computes  a 
linear  time  dependent  response  defined  by  Equations  3-1  and  3-2 
with  a,  g,  Co,  and  ko,  constants. 

White  noise  was  used  whenever  measurement  noise  was  added  to 

O 

the  signals.  The  white  noise  is  normally  distributed,  N  (0,o^) 

A  subroutine  named  NOISE  was  used  to  generate  the  noise  signal. 
The  noise  signals  were  added  to  the  generated  input  and  response 
signals  in  the  following  manner.  First,  the  excitation  and 
response  signals  were  generated  using  programs  FORCE  and  BILIN  or 
TIMVA.  Then  noise/signal  ratios  were  selected  and  used  to  obtain 
the  variances  of  the  noise  signals.  The  noise  signals  were 
generated  as  sequences  of  independent  random  variables,  and 
directly  added  to  the  excitation  and  response.  These  noisy 
signals  w^e  then  used  as  inputs  to  do  the  identification.  Note, 
no  filtering  procedure  was  used  on  the  simulated  measured  signals 
during  the  identification  process. 

Three  basic  models  to  represent  the  hysteretic  system  were 
pursued  in  this  study.  All  the  model  parameters  were  identified 
in  the  frequency  domain.  The  identification  procedures  and 
formulations  were  described  previously.  Different  identification 
approaches  may  be  applied  for  the  same  model.  Three  computer 


programs,  FREQID,  PUR3  and  PUR,  were  written  to  execute  the 
parameter  identifications. 

FREQID  performs  approximate  frequency  domain  parameter  iden¬ 
tification  for  second  and  third  order  linear  models.  It  accepts 
both  an  input  signal  from  FORCE  and  a  response  signal  from  BILIN 
or  TIMEVA.  When  desired,  the  white  noise  signals  are  added  to 
the  corresponding  input  data.  Then  FREQID  performs  the  necessary 
Fourier  transforms  and  other  data  operations.  Following  this, 
the  parameter  identification  is  executed.  One  operation 
required  in  the  parameter  identification  is  estimation  of  the 
characteristic  frequency.  This  can  be  done  simply  by  searching 
Q(a>)  +  e(u>)  for  a  minimum  value.  However,  a  more  precise  method 
for  determining  the  minimum  value  of  Q(tu)  +  e(w)  defined  in  Equa¬ 
tion  2-13  and  2-16  involves  use  of  a  least  square  method.  In 
this  improved  method,  Q(oi)  +  e(u)  terms  are  fit  into  a  polynomial 
model  where  frequency  is  the  only  variable.  Once  the  coeffi¬ 
cients  for  this  equation  are  estimated,  the  equation  representing 
the  Q(oj)  +e( (d)  is  used  to  identify  the  characteristic  frequency 
and  other  parameters. 

An  important  assumption  was  made  for  the  third-order  model, 
Equation  2-26.  In  particular,  it  was  assumed  that  a2  is  small 
in  value.  This  assumption  was  used  in  FREQID.  The  parameters  of 
this  model  may  be  identified,  without  the  assumption  that  a2  is 
small,  by  the  search  technique.  The  computer  program  PUR3  was 
written  for  this  purpose.  A  detailed  description  of  this  method 
was  given  in  Chapter  3. 

Program  PUR  performs  parameter  identification  for  second- 
order  time  varying  linear  systems.  The  approach  is  based  on  the 
procedure  described  in  Chapter  3.  The  program  accepts  the  inputs 
and  responses  generated  in  the  programs  FORCE  and  BILIN  or 
TIMEVA,  with  or  without  noise.  Four  parameters,  namely  ko,  a, 


Co.  and  g  are  identified.  The  program  employs  a  search  tech¬ 
nique;  the  initial  estimators  can  be  chosen  by  using  the  identi¬ 
fied  parameters  obtained  from  any  of  the  methods. 

Once  the  parameters  have  been  estimated,  the  energy  dissi¬ 
pated  by  the  model  is  computed.  This  result  together  with  the 
predicted  response  is  compared  to  both  the  energy  dissipated  and 
the  response  of  the  actual  system.  The  energy  and  response 
computations  are  performed  in  program  ENER2  and  ENER3  for  the 
second  and  third  order  system,  respectively. 

From  the  above  description,  the  methods  used  in  the  deter¬ 
mination  of  the  system  parameters  can  be  summarized  as  follows: 

Method  1.  Performs  parameter  identification  for  the 
second-order  system  in  the  frequency  domain 
utilizing  Equation  2-20.  No  fitting  equation  for 
Q(<d)  +  e( o>)  is  applied. 

Method  2.  Performs  parameter  identification  for  the  third- 
order  system  in  the  frequency  domain  utilizing 
the  equations  from  2-32a  to  2-35.  No  fitting 
equation  for  Q(co}  +  e(u>)  is  applied. 

Method  3.  Performs  parameter  identification  using  the  same 
approach  as  Method  1  except  the  input  data 
Q(w)  +  eu  are  replaced  by  the  fitted  polynomial 
equation.  This  additional  analysis  is  done  in  a 
subroutine  called  FIT. 

Method  4.  Performs  parameter  identification  using  the  same 
approach  as  Method  2  and  using  the  same  proce¬ 
dure  described  in  Method  3. 

Method  5.  Performs  parameter  identification  for  the  third- 
order  system  using  the  search  method  described  in 
Chapter  3.  The  operation  is  executed  in  a  pro¬ 
gram  called  PUR3.  Prior  estimates  obtained  from 
the  above  methods  are  used. 
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Method  6.  Performs  parameter  identification  for  the 

second-order  time-varying  parameter  system  in  the 
frequency  domain.  The  search  method  described  in 
Chapter  3  is  used.  This  method  also  requires 
prior  estimators  which  can  be  obtained  from  the 
information  supplied  in  Methods  1  through  5. 

All  the  methods  described  above  can  be  used  to  estimate  the 
parameters  for  the  linear  and  nonlinear  systems  even  when  noise 
is  present.  The  duration  of  the  excitation  must  be  long  enough 
to  characterize  the  system  parameters. 

In  the  following  numerical  examples,  four  basic  problems  are 
solved.  These  cases  involving  different  degrees  of  nonlinearity 
in  the  system  response  are  summarized  below. 

Case  1.  An  input  excitation  is  generated  using  Equation 

4-1.  The  input  is  used  to  excite  a  linear  SDF  sys¬ 
tem  with  viscous  damping.  The  excitation  and  linear 
response  are  used  to  identify  the  model  parameters. 
Noise  signals  can  be  added  to  the  generated  input 
and  response,  if  required. 

Case  2.  An  excitation  input  is  generated  as  in  Case  1,  but 
here  the  repsonse  of  a  bilinear  hysteretic  system  is 
computed.  Yielding  occurs  in  the  response.  The 
degree  of  nonlinearity  was  designed  using  a  compari¬ 
son  between  the  yield  displacement  of  the  bilinear 
system  and  the  maximum  displacement  of  the  linear 
system.  Let  the  yield  displacement  of  the  bilinear 
system  be  Zy.  Let  the  maximum  displacement  of  the 
linear  system  be  z^*.  In  this  case,  z^*  is 

taken  as  6.7  and  Zy  is  equal  to  6.0. 

Case  3.  Same  as  Case  2,  but  Zy  is  equal  to  5.0. 

Case  4.  Same  as  Case  2,  but  zv  is  equal  to  4.0. 


Case  5.  An  input  excitation  is  generated  as  in  Case  1.  The 
input  is  used  to  excite  a  linear  SDF  system  with 
time  varying  damping  and  stiffness.  The  excitation 
and  response  are  used  to  identify  the  model  param¬ 
eters.  Noise  signals  are  added  to  the  simulated 
input  and  response  when  required. 

4.1  Example  1 

This  example  carries  out  the  parameter  identification  using 
the  methods  described  above.  Specifically,  methods  1  through  6 
are  used  to  identify  the  parameters.  The  parameters  of  the  input 
excitation  are  listed  in  Table  4.1.  The  notation  for  the  param¬ 
eters  was  specified  above. 

TABLE  4.1.  PARAMETERS  OF  THE  FORCING  FUNCTION 

a  =  0.1  N  =  50  c j  =  10.0  j  =  1, . .  ,50 
toj  =  (1.8  +  0.008j)ir  ,j  =  1, . .  .50 
At  *  0.05  n  «  1024 

A  typical  forcing  function  history  generated  by  using  these 
parameters  is  shown  in  Figure  4.1.  Actual  forcing  functions 
measured  in  the  field  usually  contain  a  certain  amount  of  noise. 
Inputs  with  noise  to  signal  ratios  of  six  and  eight  percent  are 
shown  in  Figures  4.2  and  4.3,  respectively. 

The  response  of  some  SDF  systems  to  the  forcing  input  were 
computed.  The  energy  dissipated  in  each  structure  is  listed 
with  the  structural  parameters  in  Tables  4.1A  and  4. IB.  All 
cases  were  described  above.  The  notation  of  the  system  param¬ 
eters  is  as  follows:  k  is  the  initial  stiffness;  c  is  viscous 
damping;  ky  is  the  yield  stiffness;  Zy  is  the  yield  displace¬ 
ment;  Zfnax  is  the  maximum  displacement  of  an  SDF  system;  m  is 
the  mass  of  the  SDF  structure. 


Force  -  Kips 


Force  -  Kips 


TABLE  4.1A.  SYSTEM  PARAMETERS 


k  =  39.48  c  = 

1.257  m  = 

1.0 

Zmax  =  6* 

7  ky  =  0.0 

At  =  0.05 

n  =  1024 

TABLE  4. IB 

.  ENERGY  DISSIPATION  FOR  CASE  1 

THROUGH  CASE  4 

Cases  zv 

Damping  Energy 

Spring  Energy 

Total  Energy 

<7 

Dissipated 

Dissipated 

Dissipated 

Case  1  » 

11028.20 

0.0 

11028.2 

Case  2  6 

10199.30 

405.21 

10604.0 

Case  3  5 

8364.27 

1186.74 

3551.0 

Case  4  4 

6300.53 

1924.48 

8225.0 

First,  the  response  of  a  linear  system  (case  1)  was  com¬ 
puted.  Then,  a  slightly  nonlinear  response  of  an  SDF  structure 
was  computed  for  analysis  in  case  2.  The  displacement  response 
versus  time  for  case  2  is  plotted  in  Figure  4.4.  The  spring 
restoring  force  versus  displacement  is  shown  in  Figure  4.5.  A 
small  plastic  deformation  is  shown.  The  total  restoring  force 
versus  displacement  is  plotted  in  Figure  4.6.  The  measured 
responses  for  Case  2  with  a  certain  amount  of  noise  are  plotted 
in  Figures  4.7  and  4.8.  Figure  4.7  shows  a  measured  signal  with 
six  percent  noise  to  signal  ratio.  Figure  4.8  shows  a  measured 
signal  with  ten  percent  noise  to  signal  ratio.  Two  more  severe 
nonlinear  responses  were  computed  for  analysis  in  cases  3  and  4. 

The  displacement  response  versus  time  for  case  4  is  plotted 
in  Figure  4.9,  and  the  measured  response  for  case  4  with  ten  per¬ 
cent  noise  to  signal  ratio  is  shown  in  Figure  4.10.  The  spring 
restoring  force  versus  displacement  is  shown  in  Figure  4.11.  A 
considerable  permanent  set  is  evident  in  Figure  4.9.  Figure  4.11 
shows  that  plastic  deformation  occurs  in  the  structure  in  both 
directions  of  motion. 
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Figure  4.4  Displacement  response  of  nonlinear 
SDF  system  to  force  in  Figure  4.1 
(Case  2) 
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Using  the  forcing  function  input,  described  above,  and  the 
computed  responses,  the  parameters  of  the  structures  were  identi¬ 
fied.  The  results  of  the  parameter  identification  are  given  in 
Tables  4.2  through  4.5.  These  results  provide  the  identified 
parameters  in  the  noise-free  case.  The  energy  dissipated  by  the 
identified  system  is  listed  next  to  the  identified  parameters. 

For  method  6,  the  parameters  ao,  ai»  a2,  and  03  identify  with  the 
parameters  k0,  $,  Co,  and  a,  respectively. 

TABLE  4.2.  IDENTIFIED  PARAMETERS  AND  ENERGY 
DISSIPATED  FOR  CASE  1. 


Method 

ao 

al 

a2 

a3 

Energy 

1 

39.17 

1.26 

11210.0 

2 

39.17 

1.28 

0.0 

11090.0 

3 

40.01 

1.28 

10720.0 

4 

40.01 

2.48 

0.0236 

9270.0 

5 

39.50 

1.26 

0.0 

11200.0 

6 

39.33 

0.0 

1.274 

0.0 

10741.0 

TABLE 

4.3. 

IDENTIFIED 

PARAMETERS 

AND  ENERGY 

DISSIPATED  FOR  CASE  2 

Method 

a0 

al 

a2 

a3 

Energy 

1 

33.27 

1.45 

7968.0 

2 

33.27 

3.98 

0.05 

7028.0 

3 

35.84 

1.28 

10500.0 

4 

35.84 

3.57 

0.042 

7690.0 

5 

34.27 

3.68 

0.062 

10370.0 

6 

36.56 

0.001 

1.094 

0.022 

10077.0 

TABLE  4.4.  IDENTIFIED  PARAMETERS  AND  ENERGY 


DISSIPATED  FOR  CASE  3 


Method 

a0 

ai 

a2 

a3 

Energy 

1 

33.27 

1.35 

7802.0 

2 

33.27 

2.98 

0.032 

7368.0 

3 

35.32 

1.41 

9482.0 

4 

35.32 

2.74 

0.027 

8185.0 

5 

32.54 

2.98 

0.042 

8324.0 

6 

40.21 

-  0.0013 

1.589  - 

0.018 

9661.0 

TABLE  4.5.  IDENTIFIED  PARAMETERS  AND  ENERGY 
DISSIPATED  FOR  CASE  4 


Method 

ao 

ai 

a2 

a3 

Energy 

1 

33.27 

1.67 

7376.0 

2 

33.27 

2.28 

0.014 

7053.0 

3 

34.61 

1.75 

7703.0 

4 

34.61 

2.50 

0.018 

7486.0 

5 

30.30 

3.95 

0.066 

7598.0 

6 

40.26 

-  0.0018 

1.708 

0.001 

8186.5 

Section  4.0  demonstrated  that  when  the  higher  order  linear 
model  is  used  to  simulate  the  actual  system  behavior,  the  param¬ 
eter  a0  must  be  estimated  first  in  the  identification  procedures 
methods  1  through  4.  The  estimation  of  this  parameter  can  be 
executed  either  by  simply  searching  for  a  minimum  in  Q(w)  +  e(co) 
or  by  using  a  curve-fit  to  Q(w)  +  e(a>),  and  then  finding  the 
minimum  of  the  curve.  Figure  4.12  shows  a  realization  of  Q(u>) 
for  a  specific  case.  This  is  the  ratio  of  the  Fourier  transform 
moduli  of  the  structure  input  and  response.  An  example  of  the 
quantity  Q(o»)  +  e(w)  is  shown  in  Figure  4.13.  It 


QU) 


Frequency(w)  -  Second/cycle 
Figure  4.12  A  realization  of  Q(u) 


Frequency (u)  -  Second/cycle 
Figure  4.13  A  realization  of  Q(w)  +  e(u) 


is  apparent  from  this  diagram  why  the  use  of  a  curve-fit  provides 
better  results. 

Figures  4.14  through  4.16  show  comparisons  between  the 
responses  of  the  identified  systems  computed  by  different  meth¬ 
ods,  and  the  actual  response  of  the  bilinear  hysteretic  system  in 
case  2.  The- responses  of  the  identified  systems  match  the 
response  of  the  actual  system  so  closely  that  it  is  difficult  to 
distinguish  the  two  responses  in  these  figures.  More  identified 
responses  for  case  4  are  shown  in  Figures  4.17  through  4.19.  The 
model  responses  do  not  match  the  actual  response  as  closely  when 
residual  deformation  exists  in  the  actual  structure  since  the 
models  cannot  accumulate  permanent  deformation.  However,  peak 
responses  in  the  models  match  the  actual  system  response  quite 
well. 

In  this  example,  noise  was  added  to  the  forcing  function  and 


response  signals 

;  then  the  system  parameters 

were 

identified 

The  results  are 

summarized 

in  Tables 

4.6  and 

4.7. 

TABLE  4.6. 

IDENTIFIED 

PARAMETERS  AND  ENERGY 

DISSIPATED 

FOR  CASE 

1  WITH  TEN 

PERCENT 

NOISE  TO 

SIGNAL  RATIO 

Method 

a0 

al 

a2 

a3 

Energy 

1 

39.17 

1.37 

11022.0 

2 

39.17 

1.45 

0.003 

10900.0 

3 

39.93 

1.42 

10300.0 

4 

39.93 

2.44 

0.022 

9497.0 

5 

34.30 

2.98 

0.051 

11710.0 

6 

38.44 

0.0 

1.286 

0.0 

10929.0 

TABLE  4.7.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 


FOR  CASE 

4  WITH  TEN 

PERCENT 

NOISE  TO 

SIGNAL 

RATIO 

Method 

a0 

ai 

a2 

a3 

Energy 

1 

33.27 

1.89 

7092.0 

2 

33.27 

1.42 

0.021 

9980.0 

3 

32.02 

1.95 

6476.0 

4 

32.02 

2.60 

0.021 

6803.0 

5 

34.46 

2.70 

0.043 

11400.0 

6 

32.83 

0.0026 

1.543 

0.011 

7016.0 

These  results  show  that  the  parameter  identification  procedure  is 
still  effective  when  noise  is  present. 

4.2  Example  2 

In  this  example,  a  parameter  identification  problem  is 
solved  using  the  frequency  domain  approach.  The  methods  used  to 
identify  the  parameters  were  described  in  section  4.0,  namely 
methods  1  through  6.  The  same  forcing  function  as  illustrated  in 
Example  1  is  used.  The  only  difference  in  this  example  is  that 
the  response  was  simulated  by  a  second-order  time  varying  param¬ 
eter  system.  Unlike  the  responses  simulated  in  Example  1,  this 
example  is  a  linear  system  with  time  dependent  stiffness  and 
damping.  The  parameters  of  the  system  and  total  energy  dissi¬ 
pated  are  listed  in  Table  4.8. 

TABLE  4.8.  SYSTEM  PARAMETERS  FOR  CASE  5 

k0  *  39.48  c0  =  1.257 

8  =  -  0.01  a  =  0.01 

Total  Energy  =  9799.0 

The  definitions  of  the  symbols  are  the  same  as  in  Equation  3.2. 
This  case  was  described  in  section  4.0  as  case  5. 


The  displacement  response  versus  time  for  case  5  is  shown  in 
Figure  4.20.  The  total  restoring  force  versus  displacement  for 
this  case  is  illustrated  in  Figure  4.21.  Note  that  the  major 
axis  of  the  loops  depicted  in  the  diagram  have  different  slopes. 
This  occurs  because  the  system  stiffness  diminishes  with  time. 

The  parameters  identified  using  methods  1  through  6  together  with 
the  total  energy  dissipated  in  the  corresponding  systems  are 
shown  in  Table  4.9. 

TABLE  4.9.  IDENTIFIED  PARAMETERS  AND  ENERGY 
DISSIPATED  FOR  CASE  5 


Method 

a0 

ai 

a2 

a3 

Energy 

1 

34.70 

1.36 

9266.0 

2 

34.70 

3.25 

0.039 

7883.0 

3 

37.15 

1.29 

10960.0 

4 

37.15 

2.82 

0.031 

9074.0 

5 

34.31 

2.97 

0.05 

11160.0 

6 

37.56 

-  0.006 

1.23 

0.019 

9120.2 

shown  that 

methods 

1  through 

6  can  also 

be  used  to  identify 

model  parameters  when  noise  is  present.  The  results  obtained 
when  the  measured  signals  contain  noise  are  shown  in  Table  4.10. 

TABLE  4.10.  IDENTIFIED  PARAMETERS  AND  ENERGY  DISSIPATED 


FOR  CASE 

5  WITH 

SIX  PERCENT 

NOISE  TO 

SIGNAL 

RATIO 

Method 

a0 

al 

a2 

a3 

Energy 

1 

34.70 

1.38 

9419.0 

2 

34.70 

2.75 

0.03 

8623.0 

3 

37.34 

1.31 

11180.0 

4 

37.34 

2.24 

0.021 

10280.0 

5 

34.32 

2.49 

0.038 

11040.0 

6 

38.30 

-  0.007 

1.254 

0.025 

9245.6 
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Figure  4.20  Displacement  response  of  second-order 
time  varying  parameter  system  to  force 
in  Figure  4.1  (Case  5) 


Displacement  -  in. 

Figure  4.21  Total  restoring  force  versus  displace¬ 
ment  for  time  varying  parameter  system 


Figures  4.22,  4.24,  and  4.25  compare  the  identified  system 
responses  to  the  response  of  the  actual  system.  The  simulated 
and  actual  responses  match  quite  well  in  all  cases.  The  model 
including  time  parameters  provides  the  best  match.  Figure  4.23 
shows  the  total  restoring  force  versus  displacement.  A  change  in 
slope  of  the  major  axis  of  the  loops  is  observed.  This  behavior 
matches  the  real  system  behavior  shown  in  Figure  4.21. 
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Figure  4.24  Comparison  between  measured  response 
(heavy)  and  second-order  identified 
response  (light)  using  method  3  for 
Case  5 
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Figure  4.25  Comparison  between  measured  response 
(heavy)  and  the  third-order  response 
(light)  using  method  4  for  Case  5 
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CHAPTER  5 

SUMMARY  AND  CONCLUSIONS 

The  objective  of  this  study  was  to  develop  approximate 
linear  models  for  the  simulation  of  inelastic  system  response  and 
the  measurement  of  damage  accunul ation  in  a  structure.  It  was 
assumed  that  energy  dissipated  is  related  to  the  accumulation  of 
damage.  The  model  parameters  were  identified;  then  the  energy 
dissipated  during  a  strong  motion  was  calculated.  The  displace¬ 
ment  response  and  the  energy  dissipated  in  each  model  were  com¬ 
pared  with  the  displacement  response  and  energy  dissipated  in  the 
actual  structure. 

Three  basic  models  were  considered  in  this  study.  These  are 
second  and  third  order  linear  models  with  constant  coefficients, 
and  a  second  order  linear  model  with  time-varying  parameters. 

The  frequency  domain  approach  was  used  in  all  the  parameter  iden¬ 
tification  computations. 

Several  numerical  examples  were  solved.  The  results 
obtained  lead  to  the  following  conclusions: 

1.  Linear  and  nonlinear  hysteretic  SDF  systems  can,  in  some 
respects,  be  accurately  modeled  using  second-  and 
third-order  linear  differential  equations  with  constant 
coefficients,  and  a  second-order  linear  differential 
equation  with  time-varying  coefficients.  Specifically, 
the  models  provide  accurate  simulations  when  displace¬ 
ment  response  and  energy  dissipated  criteria  are  used. 

2.  The  frequency  domain  approach  can  be  used  to  identify 
model  parameters  of  all  three  models  when  the  force  and 
response  measuements  are  noisy. 

3.  The  second  order  model  with  time-varying  coefficients 
provides  the  best  simulation  of  system  response  and 
energy  dissipated  among  the  three  models  considered. 


While  the  procedures  developed  in  this  investigation  provide 
means  for  the  simulation  of  response  and  the  estimation  of  damage 
in  inelastic  structures,  some  improvements  can  be  made.  The 
systems  considered  in  this  study  are  SDF;  future  investigations 
should  include  multi-degree-of-freedom  structures.  The  models 
used  in  this  study  do  not  permit  the  accumulation  of  plastic 
deformation;  future  investigations  should  consider  models  that 
allow  plastic  deformation  to  accumulate.  The  tests  that  are 
summarized  in  the  appendix  show  that  material  damage  is  related 
to  energy  dissipated;  further  experiments  should  be  performed, 
and  a  mathematical  model  characterizing  the  results  should  be 
developed.  Finally,  analyses  should  be  performed  to  establish 
the  special  distribution  of  energy  dissipated  in  actual 
structural  members. 
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APPENDIX  A 

ENERGY  DISSIPATED  RELATED  TO  CONCRETE  DAMAGE 
Introduction 

This  presentation  describes  and  evaluates  an  experimental 
study  of  the  strength  reduction  and  behavior  of  plain  concrete 
subjected  to  cyclic  loading.  It  is  recognized  that  concrete  is 
damaged  by  application  of  stresses  lower  than  the  ultimate 
stress.  The  concrete  fracture  process  begins  at  very  low  stress 
and  is  continuous. 

The  damage  caused  by  loading  to  small  stresses  is  slight  and 
each  subsequent  loading  over  the  same  stress  range  produces  a 
negligible  increase  in  damage.  However,  as  the  loading  stress  is 
increased,  more  damage  occurs.  Stresses  with  peak  values  in  the 
range  of  40  percent  to  100  percent  of  the  ultimate  stress  produce 
considerable  damage  and  subsequent  loading  over  the  same  range 
cannot  be  neglected.  In  practical  situations,  when  a  severe 
excitation  is  applied  to  a  structure,  it  is  not  uncommon  for  the 
peak  stress  to  go  beyond  the  50  percent  level  of  ultimate 
stress. 

When  loading  is  repeated,  damage  accumulates  in  a  concrete 
specimen;  consequently  it  no  longer  retains  its  original 
strength.  This  concept  suggests  that  it  might  be  useful  to 
attempt  a  quantitative  evaluation  of  damage  occurring  in  concrete 
during  the  cyclic  loading.  The  objective  of  this  study  is  to 
demonstrate  that  concrete  damage  and  strength  reduction  are 
related  to  energy  dissipation  under  repeated  loading. 

When  energy  is  dissipated  during  loading  and  unloading,  an 
hysteresis  loop  is  formed  in  the  stress-strain  curve  as  shown  in 
Figure  Al.  Hie  area  enclosed  represents  the  total  energy  dissi¬ 
pated  during  one  cycle  of  loading.  This  dissipated  energy  may  be 
classified  into  two  parts,  namely,  damage  and  damping  energy 


Figure  AT.  Idealized  first  and  second  cycles 
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dissipation.  However,  the  total  accumulated  energy  dissipation 
is  of  primary  interest  in  this  experimental  study.  A  more 
detailed  discussion  of  the  energy  dissipation  mechanism  is  given 
in  Reference  [14]. 

There  are  many  methods  that  can  be  used  to  detect  and  assess 
damage  In  concrete.  Among  the  most  frequently  used  techniques 
are  those  which  assess  change  in  initial  elastic  modulus,  and 
those  which  measure  acoustic  emissions,  change  in  pulse  velocity, 
and  energy  dissipation.  In  the  present  study,  the  dissipated 
energy  method  will  be  adopted.  Attention  will  be  focused  on  this 
means  for  measuring  damage  because  of  its  intuitive  relationship 
with  the  energy  dissipation  of  an  hysteretic  structure  under 
dynamic  loads.  Other  methods  may  be  considered  as  in  References 
[14,  15,  and  16]. 

In  the  present  investigation  a  sequence  of  physical  experi¬ 
ments  was  performed.  In  each  experiment  a  concrete  cylinder 
(specifications  given  below)  was  loaded  in  uniaxial  compression. 
The  load  applied  to  each  cylinder  was  a  cyclic  load,  and  the 
energy  dissipated  was  calculated.  This  was  done  by  plotting  the 
stress  versus  strain  diagram  and  by  determining  the  area  enclosed 
within  the  hysteresis  loops.  Varying  amounts  of  energy  were 
dissipated  in  the  various  test  specimens,  and  upon  completion  of 
cyclic  testing  each  specimn  was  loaded  to  failure  in  order  to 
determine  its  residual  strength.  For  each  test  specimen,  dissi¬ 
pated  energy  and  residual  strength  were  recorded  and  the  relation 
between  these  quantities  was  established. 

More  details  of  the  testing  procedure  are  given  below, 
together  with  fresh  concrete  properties  observed  during  the  mix¬ 
ing.  These  may  provide  a  useful  reference  for  the  concretes  used 
in  the  test. 


The  concrete  specimens  tested  in  this  investigation  have  the 
mix  details  and  plastic  properties  of  fresh  concrete  as  shown  in 
Table  Al. 

TABLE  Al. 

Concrete  Mix  Details 

Aggregate  Ratio  of 
Cement  Coarse: Fine 
Type  of  Cement  W/C  Ratio  Ratio  Aggregate 

Type  1  A  0.53  4.8  60:40 

Plastic  Properties  of  Fresh  Concrete 


Mix  No. 

Slump 

(in) 

Air 

% 

Room 

Temperature 
(degrees  C) 

Unit 

Weight 

lb/ft* 

2  WE 

4 

3.5 

27 

145.8 

3  WE 

4  1/4 

4 

30 

148.96 

specimens 

were  all 

cast 

in  6-in  x  12- 

•in  steel  cylinder 

molds.  The  concrete  mix  proportions  were  constant  for  all  the 
specimens.  The  specimens  were  tested  at  a  consant  loading  rate 
of  1000  lb/sec  in  a  RIEHLE  compression  testing  machine.  The 
force  versus  strain  results  of  each  test  were  plotted  with  an  x-y 
electronic  recorder.  This  recorder  was  attached  by  an  electronic 
compressometer  which  is  properly  designed  for  this  specific 
purpose  as  shown  in  Figure  A2.  The  test  machine  was  properly 
calibrated  before  the  test. 

To  ensure  a  uniform  displacement,  of  the  specimens,  thin  sul¬ 
fur  caps  on  the  two  end  surfaces  of  the  specimens  were  employed 
and  were  allowed  to  harden  before  testing.  Specimens  were  cured 
in  water  in  the  curing  tank  at  25*C  for  14  days  and  28  days. 


Maximum 

Aggregate 

Size 

(in) 

3/4 


In  assessing  the  energy  dissipated  and  residual  strength, 
specimens  were  subjected  to  a  series  of  cycles  of  loading  and 
unloading.  The  specimens  were  loaded  up  to  a  value  in  the  range 
of  stress,  90  -  94  percent  of  the  ultimate  stress.  This  ensured 
that  damage  occurred  for  every  cycle  of  loading.  At  the  end  of 
each  cylce  of  loading  and  unloading,  the  testing  machine  was 
returned  to  a  rest  position,  and  reloading  was  commenced  immedi¬ 
ately.  To  ensure  that  concrete  characteristics  would  be  as 
nearly  uniform  as  possible,  all  the  tests  in  each  sequence  were 
run  in  one  day. 

Discussion  of  Results 

Numerous  physical  experiments  were  conducted  in  this  inves¬ 
tigation  and  characteristics  of  concrete  accumulating  damage  can 
be  derived  from  the  individual  tests  and  all  the  tests,  jointly. 
In  the  following  section  the  characteristics  of  individual  tests 
are  discussed  first;  then  damage  characteristics  related  to  the 
entire  test  sequence  are  discussed. 

A  typical  stress-strain  diagram  obtained  during  one  experi¬ 
ment  is  shown  in  Figure  A3.  A  nunber  of  characteristic  features 
can  be  extracted  from  this  result.  On  the  initial  cycle  the 
specimen  was  loaded  to  a  stress  near  its  ultimate  (95  to  98  per¬ 
cent).  It  can  be  seen  that  the  most  significant  change  in  behav 
ior  between  consecutive  loading  cycles  occurs  between  the  first 
and  second  cycle. 

The  first  loading  curve  shows  more  curvature  than  the  fol¬ 
lowing  reloading  curves  in  which  curvature  tends  to  diminish. 

The  reloading  curves  show  progressively  decreasing  slopes.  This 
may  be  attributed  to  structural  degradation  of  the  specimen. 

Another  measure  of  degradation  can  be  established  by  plot¬ 
ting  the  Initial  elastic  modulus  for  a  particular  cycle  versus 
energy  dissipation  prior  to  that  cycle.  This  is  shown  in  Figure 
A4.  As  the  energy  dissipated  gradually  increases,  the  initial 
elastic  modulus  diminishes. 


Strain  -  1n/1n  x  10 
Figure  A3.  Stress-strain  diagram  for  test  spe< 


The  above  discussion  was  based  on  one  typical  sample.  A 
similar  discussion  could  be  given  for  the  other  samples.  Figures 
A5,  A6,  and  A7  show  the  stress-strain  curves  for  some  other  spec¬ 
imens  tested  during  this  investigation. 

Some  general  characteristics  of  the  accumulation  of  damage 
in  concrete  specimens  can  be  derived  from  the  entire  collection 
of  results.  A  total  of  24  concrete  specimens  were  tested  in  this 
investigation. 

As  noted  earlier,  the  specimens  were  subjected  to  cyclic 
loadings  inducing  different  amounts  of  energy  dissipation  in  the 
various  cylinders.  Not  all  specimens  were  cycled  to  failure.  At 
least,  three  of  the  specimens  were  tested  for  the  determination 
of  the  ultimate  strength.  Other  specimens,  however,  were  cycled 
till  failure.  The  remainder  of  the  specimens  were  cycled  till  a 
certain  amount  of  energy  was  dissipated;  then  these  were  loaded 
to  failure  in  order  to  find  their  residual  strength. 

Using  these  data,  a  characteristic  of  the  specimens  can  be 
extracted.  The  total  energy  dissipated  by  each  particular  speci¬ 
men  was  plotted  against  the  residual  strength  of  the  specimen  as 
shown  in  Figure  A9.  Another  result  can  be  illustrated  by  plot¬ 
ting  the  total  energy  dissipated  versus  percentage  of  decrease  in 
strength  as  shown  in  Figure  A8.  Both  diagrams  show  a  decrease  in 
strength  as  the  total  energy  dissipated  is  increased.  Since  only 
a  limited  number  of  specimens  were  tested,  no  direct  mathematical 
expression  relating  the  residual  strength  to  the  total  energy 
dissipated  was  obtained.  While  such  a  relation  could  be  estab¬ 
lished,  further  testing  is  required  to  derive  a  general  relation¬ 
ship.  However,  the  pesent  results  provide  the  information  needed 
to  conclude  that  energy  dissipation  is  truly  related  to  the 
residual  strength. 
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Strain  -  in/in  x  10 

Figure  A5.  Typical  stress-strain  diagram  for  a  con¬ 
crete  test  specimen  under  cyclic  load. 
Peak  stress  90  percent  of  failure  stress 
for  Batch  #3  (14  days  curing) 


Strain  -  in/in  x  10 

Figure  A6.  Typical  stress-strain  diagram  for  a 
concrete  test  specimen  under  cyclic 
load.  Peak  stress  90  percent  of 
failure  stress  for  Batch  #2 
(28  days  curing) 


Residual  force  -  lb  x  10  Percentage  decrease  in  strength 


The  experimental  technique  described  above  provides  an 
approach  for  the  estimation  of  damage  in  concrete.  Based  on  the 
physical  experiments,  the  following  conclusions  can  be  made: 

1.  The  most  significant  change  in  the  properties  of  the 
concrete  occur  between  the  first  cycle  and  second  cycle 
when  loading  in  the  first  cycle  is  severe. 

2.  The  initial  elastic  modulus  of  the  specimens  gradually 
diminishes  as  energy  is  dissipated.  This  implies  that 
the  damage  of  concrete  under  cyclic  loading  occurs 
progressively. 

3.  The  energy  dissipated  in  a  concrete  specimen  is 
adversely  related  to  residual  strength.  As  the  energy 
dissipated  increases,  the  residual  strength  decreases. 
Therefore,  energy  dissipated  may  be  used  to  predict  the 
damage  of  a  structure  under  a  severe  loading.  Moreover 
total  energy  dissipated  may  be  considered  as  an 
indicator  of  the  degree  of  damage  in  an  hysteretic 
structure. 

Some  restrictions  apply  to  the  above  conclusions.  The  work 
is  limited  to  beha\ior  in  uniaxial  compression.  Other  types 
of  loading  are  possible  and  further  tests  are  required  to 
characterize  damage  under  general  loading.  It  has  been  assumed 
that  the  creep  effect  is  small  enough  to  be  neglected  in  this 
Investigation. 


APPENDIX  B 
COMPUTER  PROGRAM  PUR 

***************************************** 

***********  PROGRAM  NAME  "PUR''************ 
***************************************** 

THIS  PROGRAM  IDENTIFIES  THE  PARAMETERS  OF  A 
TIME  VARYING  SECOND-ORDER  LINEAR  MODEL  BY 
USING  THE  PERTURBATION  METHOD  AND  ITERATIVE 
NEWTON-  RAPHSON  PROCEDURE.  IT  MAY  ALSO  USE 
THE  POLYNOMIAL  FITTING  APPROACH. 

COMPLEX  FF( 1024) ,FD( 1024) 

COMPLEX  ZZ ( 1024 ) , HD, Z ( 1024) , HABS ( 512 ) 
DIMENSION  DD( 1024),F(1024), HMOD (512), HZ ( 512 ) 
DIMENSION  EK( 3 ) , COE ( 4 ) , S(9) , INDEX ( 9 ) , AS ( 9 ) 
COMMON/ARRAY/FF , FD 
CALL  OPSYS( 'ALLOC' ,' FD1' , 10) 

CALL  OPSYS ( ' ALLOC ' , ' RH1 ’,7) 

CALL  OPSYS ( ' ALLOC ' , ' NFL ' , 4 ) 

EXTERNAL  MATINV 

SM  :  SYSTEM  MASS 

ALPHA  :  DAMPING  COEFFICIENT 

BETA  :  STIFFNESS  COEFFICIENT 

ZETA  :  DAMPING  RATIO 

Cl  :  INITIAL  GUESS  'DAMPING' 

C2  :  INITIAL  GUESS  'ALPHA*C1' 

C3  :  INITIAL  GUESS  'STIFFNESS' 

C4  :  INITIAL  GUESS  ' BETA*C3 ' 

ESLON  :  ACCURACY  MEASURE (FOR  ESL0N1-ESL0N4 ) 
W1-W4  :  ACCURACY  MEASURE  FOR  NEWTON  METHOD 
DA1-DA4  :  INCREMENT  VALUES 

SM=1 . 0 
Cl=l .  5 
C2=0 . 02 
C3=37 . 5 
C4=-0 . 01 
ZETA=0 . 1 
ESLON1=0 . 02 
ESL0N2=0 . 02 
ESLON3=0 . 02 
ESLON4=0 . 02 
WN=SQRT ( C3/SM) 

DA1=0 . 01 
DA2=0 . 001 
DA3=1 . 0 
DA4=0 . 01 
Wl=0.01 
W2=0 . 01 
W3=0 . 01 
W4=0 . 01 

PI=3. 1415926535 
DT=0 . 05 


N=2**M 
NT=N/8 
TT=DT*N 
DW=2 . *PI/TT 
NS=N/2 
CA=P I/NS 

READ  IN  THE  TIME  DOMAIN  RESPONSE 

READ(10,80)  (F(I),DD(I), I=1,N) 

WINDOW  THE  DATA 

DO  3  1=1, NS 
IM=N-I+1 

CMULT=0 . 5* ( 1 . O-COS ( ( I-1)*CA) ) 
F(IM)=CMULT*F(IM) 

DD ( I M ) =CMULT  *  DD ( I M ) 

3  CONTINUE 

FFT  THE  INPUT  AND  RESPONSE 
DO  8  1=1, N 

FF ( I ) =CMPLX( F ( I ) , 0 . 0 ) *TT 
ZZ ( I )=CMPLX(DD( I ) , 0 . 0 ) *TT 

8  CONTINUE 

CALL  FFT1(FF,M,N, -1.0) 

CALL  FFT1 ( ZZ, M, N, -1 . 0 ) 

FIND  THE  MOD ( H ) 

DO  9  1=1, NT 
FA=CABS ( FF ( I ) ) 

ZA=CABS ( ZZ ( I ) ) 

HMOD ( I ) =ZA/FA 

9  CONTINUE 

CALCULATE  THE  DERIVATIVE  OF  FF 


DO  12  1=1, N 
TC=( 1-1 ) *DT 
F ( I )=F ( I ) *TC 

FD( I )=CMPLX( F( I ) , 0 . 0 ) *TT 

12  CONTINUE 

CALL  FFT1 ( FD, M, N, -1. ) 

DO  13  1=1, N 

13  FD ( I ) =FD ( I ) *CM?LX ( 0 . 0 , - 1 . ) 


FIND  THE  NUMBER  OF  PTS  NEEDED  IN  THE  3 AND 
NC= (1.0*  ZET A  *  WN ) /DW 


noo  noon 


NW=WN/DW 
NS=NW-NC 
NF=NW+NC 
AA1=C1 
AA2=C2 
AA3=C3 
AA4=C4 
26  L=1 

28  WRITE ( 6,89)01,02,03, C4 
IC0UNT=2 

I STEA= I COUNT 

AC1=C1 

AC2=C2 

AC3=C3 

AC4=C4 

CALL  SUBROUTINE  PURT  TO  FIND  THE  MODULUS 
OF  H  FROM  THE  MEASURED  INPUT 

29  WN=SQRT ( C3/SM ) 

NC= ( 1 . 0*ZETA*WN ) /DW 

NW=WN/DW 

NS=NW-NC 

NF=NW+NC 

CALL  PURT(DW, NT,  Cl , C2 , C3 , C4, HZ ) 

FIND  E( ERROR  MEASURE) 

EE=0 . 0 

DO  30  I=NS,NF 
E=HMOD ( I ) -HZ ( I ) 

EE=EE+(E**2) 

30  CONTINUE 

EK ( I COUNT ) =EE 
I F ( I COUNT . EQ . 3 )  GO  TO  50 
GO  TO  (31, 32,33,34) ,L 

31  IF( ICOUNT . EQ. 1 )  GO  TO  40 
C1=C1-DA1 
IC0UNT=IC0UNT-1 

GO  TO  29 

40  C1=AC1+DA1 
ICOUNT= I STEA+ 1 
GO  TO  29 

32  IF( ICOUNT. EQ.l)  GO  TO  41 
C2=C2-DA2 
ICOUNT=ICOUNT- 1 

GO  TO  29 

41  C2=AC2+DA2 
IC0UNT=ISTEA+1 
GO  TO  29 

33  IF( ICOUNT. EQ.l)  GO  TO  42 
C3=C3 -DA3 


I COUNT= I COUNT - 1 
GO  TO  29 

42  C3=AC3+DA3 
ICOUNT=ISTEA+l 
GO  TO  29 

34  IF( ICOUNT.EQ. 1)  GO  TO  43 
C4=C4-DA4 
I COUNT= I COUNT - 1 
GO  TO  29 

43  C4=AC4+DA4 
ICOUNT=ISTEA+l 
GO  TO  29 

CALCULATE  THE  FIRST  AND  SECOND  DERIVETIVES 

50  GO  T0( 60 ,  61 , 62 , 63  ) ,  L 

60  TD=DA1 
AV=AC1 
W=W1 

GO  TO  70 

61  TD=DA2 
AV=AC2 
W=W2 

GO  TO  70 

62  TD=DA3 
AV=AC3 
W=W3 

GO  TO  70 

63  TD=DA4 
AV=AC4 
W=W4 

GO  TO  70 

70  DEK=( EK( 3 ) -EK( 1 ) )/( 2 . 0*TD ) 

DDEK=(EK( 3 ) -2 . *EK( 2 ) +EK( 1) )/(TD**2) 

DO  THE  NEWTON  RAPHSON  METHOD 

ACC=AV- ( DEK/DDEK ) 

WE=ABS( (ACC-AV)/ACC) 

WD=WE-W 

IF(WD.LE.O.O)  GO  TO  87 
99  GO  TO( 101,102, 103  , 104 ), L 

101  C1=ACC 
GO  TO  28 

102  C2=ACC 
GO  TO  28 

103  C3=ACC 
GO  TO  28 

104  C4=ACC 
GO  TO  28 

87  GO  TO( 92 , 93 , 94, 95 ) , L 
92  C1=ACC 


g 
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b 
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% 
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>1%’ 
LVV-" 


GO  TO  85 

93  C2=ACC 
GO  TO  85 

94  C3=ACC 
GO  TO  85 

95  C4=ACC 
85  L=L+1 

IF(L.EQ. 5)  GO  TO  100 
GO  TO  28 
C 

C  COMPARE  THE  PARAMETERS 

C 

100  CD1=ABS ( { AA1-C1 )/Cl ) 

CD2=ABS{ ( AA2-C2 )/C2 ) 

CD3=ABS ( (AA3-C3 )/C3 ) 

CD4=ABS( ( AA4-C4)/C4) 

IF(CD1 . GT . ESL0N1 )  GO  TO  81 
IF(CD2 . GT. ESL0N2 )  GO  TO  81 
IF(CD3 . GT. ESL0N3 )  GO  TO  81 
IF( CD3 . GT . ESL0N4)  GO  TO  81 
WRITE ( 6, * )  01,02,03, C4 
GO  TO  88 

81  AA1=C1 
AA2=C2 
AA3=C3 
AA4=C4 
GO  TO  26 

80  FORMAT (2E 12. 4) 

82  FORMAT (3E13. 4) 

89  FORMAT (4E 12. 4) 

150  FORMAT (2E1 6 . 6 ) 

88  STOP 
END 
C 

C  PERTURBATION  METHOD 

C 

SUBROUTINE  PURT(DW, NT, Cl , 02 , 03 , 04, HA) 
COMPLEX  FF(1024) ,FD(1024) , HABS (512) ,2(512) 
COMPLEX  HD , H , HDEV , HI , H2 , H3 , H4 
DIMENSION  HA (NT) 

COMMON/ ARRAY/FF, FD 
SM=1 . 0 

DO  14  1=1, NT 
W=(I-1)*DW 

HD=CMPLX( (C3-SM* ( W**2 ) ) ,W*C1) 

H=1 . O/HD 

HDEV=CMPLX( (2. *SM*W) , -Cl)/ (HD* *2) 

H1=H*FF ( I ) 

H2=C2*H* (Hl+W* (HDEV*FF ( I ) +H*FD( I ) ) ) 
H3=CMPLX( 0 .0,-04) *K* (KDEV*FF ( I ) +H*FD( I ) ) 
2(1 )=H1+H2 +H3 
H4=HDEV+H* ( FD { I ) /FF ( I ) ) 


HABS ( I ) =H* ( 1 . 0+C2  * ( H+W*H4 ) -CMP LX ( 0 . 0 , C4 ) *H4 ) 
HA ( I ) =CABS ( HABS ( I ) ) 

14  CONTINUE 
GO  TO  10 
10  RETURN 
END 


FFT  PROGRAM 

******************************************** 
MULTIPLIED  BY  T  WHEN  USING  THE  FORWARD  FFT 
DIVIDED  BY  DT  FOR  BACKWARD  FFT 
SIGN  =-l.  FOR  DFT,  SIGN=1.  FOR  IDFT 
********************************************* 


SUBROUTINE  FFT1 ( A, N, NB, S IGN) 
COMPLEX  A ( NB ) ,  U ,  W ,  T 

DIVIDE  ALL  ELEMENT  BY  NB 

DO  1  J=1 , NB 

1  A(J)=A(J)/NB 

REORDER  SEQUENCE 

NBD2=NB/2 

NBM1=NB-1 

J=1 

DO  4  L=1 , NBM1 
IF(L.GE.J)  GO  TO  2 
T=A(J) 

A ( J ) —A ( L ) 

A ( L ) =T 

2  K=NBD2 

3  IF(K.GE.J)  GO  TO  4 
J=J-K 

K=K/2 
GO  TO  3 

4  J=J+K 

CALCULATE  FFT 

PI=3. 141592653589793 
DO  6  M=1 , N 
U=(l. 0,0.0) 

ME=2**M 

K=ME/2 

W=CMPLX(COS(PI/K) , SIGN*SIN(PI/K) ) 
DO  6  J=1.K 


ooo  oooooon 


5 

6 


10 


10 


A( LPK) =A( L ) -T 

A(L)=A(L)+T 

U=U*W 

RETURN 

END 

LEAST  SQUARE  POLYNOMINAL  FITTING 

A  POLYNOMIAL  FIT  SUBROUTINE 

THIS  PROGRAM  IS  PREPARED  SO  THAT  IT  CAN 

BE  USED  WHENEVER  IT  IS  NECESSARY. 

SUBROUTINE  FIT ( Y, MP , N, XI , X2 , A ) 

DIMENSION  XX(9,4),Y(9) ,X(9),XP(4,4),A(4) ,D(4,4) ,B(4) 
DIMENSION  EPS ( 9 ) 

DX= ( X2 -XI ) /FLOAT ( N- 1 ) 

DO  10  1=1  ,N 

X( I ) =X1+FL0AT ( I - 1 ) *DX 

DO  1  1=1 ,N 

XX(I<1)=1.0 

DO  1  J=2 , MP 

XX ( I , J)=X( I)**(J-1) 

DO  2  1=1 ,MP 
DO  2  J=I , MP 
XP ( I , J ) =0 . 0 
DO  2  K=1 , N 

XP( I , J)=XP( I , J)+XX(K, I ) *XX(K, J) 

DO  3  1=2 ,MP 
IM=I - 1 
DO  3  J=1 , IM 
XP(I, J)=XP(J, I) 

CALL  MATINV ( XP, MP , D ) 

DO  4  1=1, MP 
B ( I ) =0 . 0 
DO  4  J=1,N 

B( I )=B( I ) +XX ( J, I ) *Y( J ) 

DO  5  1=1, MP 

A ( I ) =0 . 0 

DO  5  J=1 , MP 

A( I )=A( I)+D(I,J)*B(J) 

RETURN 

END 

MATRIX  INVERSION 

SUBROUTINE  MATINV ( C, N, D) 

MATRIX  INVERSION  C- INPUT  D-OUTPUT 
DIMENSION  C(4,4) ,D(4,4) 

DO  10  J=1 , N 
DO  10  K=1 , N 
D( J, K)=0 . 0 
DO  11  K=1 , N 
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11  D(K, K)=l . 0 
Pl=1.0 
DO  55  1=1, N 
P2=C( 1,1) 

DO  40  J=1,N 
C(I, J)=C(I, J)/P2 
40  D ( I , J ) =D ( I , J ) /P2 
DO  51  IC=1 , N 
P3=-C(IC, I) 

DO  50  K=1,N 
IF(  IC-I ) 21, 51, 21 
21  C( IC, K)=C( I ,  K) *P3+C( IC, K) 

50  D(IC,K)=D(I,K)*P3+D(IC,K) 

51  CONTINUE 
P1=P2*P1 

IF( ( 1+2 ) -N) 55 ,53,55 
53  DET=P1*( <C(I+1, I+l)*C(I+2, 1+2) ) 
*-(C(I+2, I+1)*C(I+1, 1+2) ) ) 

55  CONTINUE 

DO  70  IT=1,N 
DO  70  IS=1,N 
70  C( IT, IS)=D( IT, IS) 

RETURN 
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